home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / bzr_sym.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  13.5 KB  |  473 lines
  1. /******************************************************************************
  2. * Bzr_Sym.c - Bezier symbolic computation.                      *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Sep. 92.                          *
  5. ******************************************************************************/
  6.  
  7. #include "cagd_loc.h"
  8.  
  9. /******************************************************************************
  10. * Given two bezier curves - multiply them coordinatewise.              *
  11. * The two curves are promoted to same point type before the multiplication    *
  12. * can take place.                                  *
  13. * Return a bezier curve representing their product.                  *
  14. ******************************************************************************/
  15. CagdCrvStruct *BzrCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
  16. {
  17.     CagdBType IsNotRational;
  18.     int i, j, k, MaxCoord, Order,
  19.     Order1 = Crv1 -> Order,
  20.     Order2 = Crv2 -> Order,
  21.     Degree1 = Order1 - 1,
  22.     Degree2 = Order2 - 1;
  23.     CagdCrvStruct *ProdCrv;
  24.     CagdRType **Points, **Points1, **Points2;
  25.  
  26.     if (!CAGD_IS_BEZIER_CRV(Crv1) || !CAGD_IS_BEZIER_CRV(Crv2))
  27.     FATAL_ERROR(CAGD_ERR_BZR_CRV_EXPECT);
  28.  
  29.     Crv1 = CagdCrvCopy(Crv1);
  30.     Crv2 = CagdCrvCopy(Crv2);
  31.     if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE))
  32.     FATAL_ERROR(CAGD_ERR_CRV_FAIL_CMPT);
  33.  
  34.     ProdCrv = BzrCrvNew(Order = Order1 + Order2 - 1, Crv1 -> PType);
  35.     IsNotRational = !CAGD_IS_RATIONAL_CRV(ProdCrv);
  36.     MaxCoord = CAGD_NUM_OF_PT_COORD(ProdCrv -> PType);
  37.  
  38.     Points = ProdCrv -> Points;
  39.     Points1 = Crv1 -> Points;
  40.     Points2 = Crv2 -> Points;
  41.  
  42.     for (i = 0; i < Order; i++)
  43.     for (k = IsNotRational; k <= MaxCoord; k++)
  44.         Points[k][i] = 0.0;
  45.  
  46.     for (i = 0; i < Order1; i++) {
  47.     for (j = 0; j < Order2; j++) {
  48.         CagdRType
  49.         Coef = CagdIChooseK(i, Degree1) *
  50.                CagdIChooseK(j, Degree2) /
  51.                CagdIChooseK(i + j, Degree1 + Degree2);
  52.  
  53.         for (k = IsNotRational; k <= MaxCoord; k++)
  54.         Points[k][i+j] += Coef * Points1[k][i] * Points2[k][j];
  55.     }
  56.     }
  57.  
  58.     CagdCrvFree(Crv1);
  59.     CagdCrvFree(Crv2);
  60.  
  61.     return ProdCrv;
  62. }
  63.  
  64. /******************************************************************************
  65. * Given two bezier curve lists - multiply them one at a time.              *
  66. * Return a bezier curve lists representing their products.              *
  67. ******************************************************************************/
  68. CagdCrvStruct *BzrCrvMultList(CagdCrvStruct *Crv1Lst, CagdCrvStruct *Crv2Lst)
  69. {
  70.     CagdCrvStruct *Crv1, *Crv2,
  71.     *ProdCrvTail = NULL,
  72.     *ProdCrvList = NULL;
  73.  
  74.     for (Crv1 = Crv1Lst, Crv2 = Crv2Lst;
  75.      Crv1 != NULL && Crv2 != NULL;
  76.      Crv1 = Crv1 -> Pnext, Crv2 = Crv2 -> Pnext) {
  77.     CagdCrvStruct
  78.         *ProdCrv = BzrCrvMult(Crv1, Crv2);
  79.  
  80.     if (ProdCrvList == NULL)
  81.         ProdCrvList = ProdCrvTail = ProdCrv;
  82.     else {
  83.         ProdCrvTail -> Pnext = ProdCrv;
  84.         ProdCrvTail = ProdCrv;
  85.     }
  86.     }
  87.  
  88.     return ProdCrvList;
  89. }
  90.  
  91. /******************************************************************************
  92. * Given two bezier surfaces - multiply them coordinatewise.              *
  93. * The two surfaces are promoted to same point type before the multiplication  *
  94. * can take place.                                  *
  95. * Return a bezier surface representing their product.                  *
  96. ******************************************************************************/
  97. CagdSrfStruct *BzrSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
  98. {
  99.     CagdBType IsNotRational;
  100.     int i, j, k, l, m, MaxCoord, UOrder, VOrder,
  101.     UOrder1 = Srf1 -> UOrder,
  102.     VOrder1 = Srf1 -> VOrder,
  103.     UOrder2 = Srf2 -> UOrder,
  104.     VOrder2 = Srf2 -> VOrder,
  105.     UDegree1 = UOrder1 - 1,
  106.     VDegree1 = VOrder1 - 1,
  107.     UDegree2 = UOrder2 - 1,
  108.     VDegree2 = VOrder2 - 1;
  109.     CagdSrfStruct *ProdSrf;
  110.     CagdRType **Points, **Points1, **Points2;
  111.  
  112.     if (!CAGD_IS_BEZIER_SRF(Srf1) || !CAGD_IS_BEZIER_SRF(Srf2))
  113.     FATAL_ERROR(CAGD_ERR_BZR_SRF_EXPECT);
  114.  
  115.     Srf1 = CagdSrfCopy(Srf1);
  116.     Srf2 = CagdSrfCopy(Srf2);
  117.     if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE))
  118.     FATAL_ERROR(CAGD_ERR_SRF_FAIL_CMPT);
  119.  
  120.     ProdSrf = BzrSrfNew(UOrder = UOrder1 + UOrder2 - 1,
  121.             VOrder = VOrder1 + VOrder2 - 1, Srf1 -> PType);
  122.     IsNotRational = !CAGD_IS_RATIONAL_SRF(ProdSrf);
  123.     MaxCoord = CAGD_NUM_OF_PT_COORD(ProdSrf -> PType);
  124.  
  125.     Points = ProdSrf -> Points;
  126.     Points1 = Srf1 -> Points;
  127.     Points2 = Srf2 -> Points;
  128.  
  129.     for (k = IsNotRational; k <= MaxCoord; k++) {
  130.     CagdRType
  131.         *Pts = Points[k];
  132.  
  133.     for (i = 0; i < UOrder * VOrder; i++)
  134.         *Pts++ = 0.0;
  135.     }
  136.  
  137.     for (i = 0; i < UOrder1; i++) {
  138.     for (j = 0; j < VOrder1; j++) {
  139.         for (l = 0; l < UOrder2; l++) {
  140.         for (m = 0; m < VOrder2; m++) {
  141.             int Index = CAGD_MESH_UV(ProdSrf, i + l, j + m),
  142.             Index1 = CAGD_MESH_UV(Srf1, i, j),
  143.             Index2 = CAGD_MESH_UV(Srf2, l, m);
  144.             CagdRType
  145.             Coef1 = CagdIChooseK(i, UDegree1) *
  146.                     CagdIChooseK(l, UDegree2) /
  147.                 CagdIChooseK(i + l, UDegree1 + UDegree2),
  148.             Coef2 = CagdIChooseK(j, VDegree1) *
  149.                     CagdIChooseK(m, VDegree2) /
  150.                 CagdIChooseK(j + m, VDegree1 + VDegree2);
  151.  
  152.             for (k = IsNotRational; k <= MaxCoord; k++)
  153.             Points[k][Index] += Coef1 * Coef2 * 
  154.                 Points1[k][Index1] * Points2[k][Index2];
  155.         }
  156.         }
  157.     }
  158.     }
  159.  
  160.     CagdSrfFree(Srf1);
  161.     CagdSrfFree(Srf2);
  162.  
  163.     return ProdSrf;
  164. }
  165.  
  166. /******************************************************************************
  167. * Given a rational bezier curve - computes its derivative curve using the     *
  168. * quotient rule for differentiation.                          *
  169. ******************************************************************************/
  170. CagdCrvStruct *BzrCrvDeriveRational(CagdCrvStruct *Crv)
  171. {
  172.     CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *DCrvW, *DCrvX, *DCrvY, *DCrvZ,
  173.     *TCrv1, *TCrv2, *DeriveCrv;
  174.  
  175.     CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
  176.     if (CrvW)
  177.     DCrvW = BzrCrvDerive(CrvW);
  178.     else {
  179.         DCrvW = NULL;
  180.     FATAL_ERROR(CAGD_ERR_RATIONAL_EXPECTED);
  181.     }
  182.     
  183.     if (CrvX) {
  184.     DCrvX = BzrCrvDerive(CrvX);
  185.  
  186.     TCrv1 = BzrCrvMult(DCrvX, CrvW);
  187.     TCrv2 = BzrCrvMult(CrvX, DCrvW);
  188.  
  189.     CagdCrvFree(CrvX);
  190.     CagdCrvFree(DCrvX);
  191.     CrvX = CagdCrvSub(TCrv1, TCrv2);
  192.     CagdCrvFree(TCrv1);
  193.     CagdCrvFree(TCrv2);
  194.     }
  195.  
  196.     if (CrvY) {
  197.     DCrvY = BzrCrvDerive(CrvY);
  198.  
  199.     TCrv1 = BzrCrvMult(DCrvY, CrvW);
  200.     TCrv2 = BzrCrvMult(CrvY, DCrvW);
  201.  
  202.     CagdCrvFree(CrvY);
  203.     CagdCrvFree(DCrvY);
  204.     CrvY = CagdCrvSub(TCrv1, TCrv2);
  205.     CagdCrvFree(TCrv1);
  206.     CagdCrvFree(TCrv2);
  207.     }
  208.  
  209.     if (CrvZ) {
  210.     DCrvZ = BzrCrvDerive(CrvZ);
  211.  
  212.     TCrv1 = BzrCrvMult(DCrvZ, CrvW);
  213.     TCrv2 = BzrCrvMult(CrvZ, DCrvW);
  214.  
  215.     CagdCrvFree(CrvZ);
  216.     CagdCrvFree(DCrvZ);
  217.     CrvZ = CagdCrvSub(TCrv1, TCrv2);
  218.     CagdCrvFree(TCrv1);
  219.     CagdCrvFree(TCrv2);
  220.     }
  221.  
  222.     TCrv1 = BzrCrvMult(CrvW, CrvW);
  223.     CagdCrvFree(CrvW);
  224.     CrvW = TCrv1;
  225.  
  226.     if (!CagdMakeCrvsCompatible(&CrvW, &CrvX, TRUE, TRUE) ||
  227.     !CagdMakeCrvsCompatible(&CrvW, &CrvY, TRUE, TRUE) ||
  228.     !CagdMakeCrvsCompatible(&CrvW, &CrvZ, TRUE, TRUE))
  229.     FATAL_ERROR(CAGD_ERR_CRV_FAIL_CMPT);
  230.  
  231.     DeriveCrv = CagdCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ);
  232.  
  233.     if (CrvX)
  234.     CagdCrvFree(CrvX);
  235.     if (CrvY)
  236.     CagdCrvFree(CrvY);
  237.     if (CrvZ)
  238.     CagdCrvFree(CrvZ);
  239.     if (CrvW) {
  240.     CagdCrvFree(CrvW);
  241.     CagdCrvFree(DCrvW);
  242.     }
  243.  
  244.     return DeriveCrv;
  245. }
  246.  
  247. /******************************************************************************
  248. * Given a rational bezier surface - computes its derivative curve using the   *
  249. * quotient rule for differentiation.                          *
  250. ******************************************************************************/
  251. CagdSrfStruct *BzrSrfDeriveRational(CagdSrfStruct *Srf, CagdSrfDirType Dir)
  252. {
  253.     CagdSrfStruct *SrfW, *SrfX, *SrfY, *SrfZ, *DSrfW, *DSrfX, *DSrfY, *DSrfZ,
  254.     *TSrf1, *TSrf2, *DeriveSrf;
  255.  
  256.     CagdSrfSplitScalar(Srf, &SrfW, &SrfX, &SrfY, &SrfZ);
  257.     if (SrfW)
  258.     DSrfW = BzrSrfDerive(SrfW, Dir);
  259.     else {
  260.         DSrfW = NULL;
  261.     FATAL_ERROR(CAGD_ERR_RATIONAL_EXPECTED);
  262.     }
  263.  
  264.     if (SrfX) {
  265.     DSrfX = BzrSrfDerive(SrfX, Dir);
  266.  
  267.     TSrf1 = BzrSrfMult(DSrfX, SrfW);
  268.     TSrf2 = BzrSrfMult(SrfX, DSrfW);
  269.  
  270.     CagdSrfFree(SrfX);
  271.     CagdSrfFree(DSrfX);
  272.     SrfX = CagdSrfSub(TSrf1, TSrf2);
  273.     CagdSrfFree(TSrf1);
  274.     CagdSrfFree(TSrf2);
  275.     }
  276.     if (SrfY) {
  277.     DSrfY = BzrSrfDerive(SrfY, Dir);
  278.  
  279.     TSrf1 = BzrSrfMult(DSrfY, SrfW);
  280.     TSrf2 = BzrSrfMult(SrfY, DSrfW);
  281.  
  282.     CagdSrfFree(SrfY);
  283.     CagdSrfFree(DSrfY);
  284.     SrfY = CagdSrfSub(TSrf1, TSrf2);
  285.     CagdSrfFree(TSrf1);
  286.     CagdSrfFree(TSrf2);
  287.     }
  288.     if (SrfZ) {
  289.     DSrfZ = BzrSrfDerive(SrfZ, Dir);
  290.  
  291.     TSrf1 = BzrSrfMult(DSrfZ, SrfW);
  292.     TSrf2 = BzrSrfMult(SrfZ, DSrfW);
  293.  
  294.     CagdSrfFree(SrfZ);
  295.     CagdSrfFree(DSrfZ);
  296.     SrfZ = CagdSrfSub(TSrf1, TSrf2);
  297.     CagdSrfFree(TSrf1);
  298.     CagdSrfFree(TSrf2);
  299.     }
  300.     
  301.     TSrf1 = BzrSrfMult(SrfW, SrfW);
  302.     CagdSrfFree(SrfW);
  303.     SrfW = TSrf1;
  304.  
  305.     if (!CagdMakeSrfsCompatible(&SrfW, &SrfX, TRUE, TRUE, TRUE, TRUE) ||
  306.     !CagdMakeSrfsCompatible(&SrfW, &SrfY, TRUE, TRUE, TRUE, TRUE) ||
  307.     !CagdMakeSrfsCompatible(&SrfW, &SrfZ, TRUE, TRUE, TRUE, TRUE))
  308.     FATAL_ERROR(CAGD_ERR_SRF_FAIL_CMPT);
  309.  
  310.     DeriveSrf = CagdSrfMergeScalar(SrfW, SrfX, SrfY, SrfZ);
  311.  
  312.     if (SrfX)
  313.     CagdSrfFree(SrfX);
  314.     if (SrfY)
  315.     CagdSrfFree(SrfY);
  316.     if (SrfZ)
  317.     CagdSrfFree(SrfZ);
  318.     if (SrfW) {
  319.     CagdSrfFree(SrfW);
  320.     CagdSrfFree(DSrfW);
  321.     }
  322.  
  323.     return DeriveSrf;
  324. }
  325.  
  326. /******************************************************************************
  327. * Given a bezier curve - convert it to (possibly) piecewise cubic.          *
  328. * If the curve is                                  *
  329. * 1. A cubic - a copy if it is returned.                      *
  330. * 2. Lower than cubic - a degree raised (to a cubic) curve is returned.       *
  331. * 3. Higher than cubic - a C^1 piecewise cubic approximation is computed.     *
  332. *    Note in this case a list of polynomial cubic curves is returned. Tol is  *
  333. *    then used for the distance tolerance error measure for the approximation *
  334. *   If, however, NoRational is set, rational curves of any order will also be *
  335. * approximated using cubic polynomials. Furthermore if the total length of    *
  336. * control polygon is more than MaxLen, the curve is subdivided until this is  *
  337. * not the case.                                      *
  338. ******************************************************************************/
  339. CagdCrvStruct *BzrApproxBzrCrvAsCubics(CagdCrvStruct *Crv, CagdRType Tol,
  340.                        CagdRType MaxLen, CagdBType NoRational)
  341. {
  342.     CagdCrvStruct *TCrv, *CubicCrvs,
  343.     *CubicCrvsMaxLen = NULL;
  344.  
  345.     if (CAGD_IS_RATIONAL_CRV(Crv) && NoRational)
  346.     return BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol);
  347.  
  348.     switch (Crv -> Order) {
  349.     case 2:
  350.         CubicCrvs = BzrCrvDegreeRaiseN(Crv, 4);
  351.         break;
  352.     case 3:
  353.         CubicCrvs = BzrCrvDegreeRaise(Crv);
  354.         break;
  355.     case 4:
  356.         CubicCrvs = CagdCrvCopy(Crv);
  357.     default:
  358.         if (Crv -> Order < 4)
  359.         FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
  360.         CubicCrvs = BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol);
  361.     }
  362.  
  363.     for (TCrv = CubicCrvs; TCrv != NULL; TCrv = TCrv -> Pnext) {
  364.     CagdCrvStruct *TCrv2,
  365.         *MaxLenCrv = CagdLimitCrvArcLen(TCrv, MaxLen);
  366.  
  367.     for (TCrv2 = MaxLenCrv;
  368.          TCrv2 -> Pnext != NULL;
  369.          TCrv2 = TCrv2 -> Pnext);
  370.  
  371.     TCrv2 -> Pnext = CubicCrvsMaxLen;
  372.     CubicCrvsMaxLen = MaxLenCrv;
  373.     }
  374.  
  375.     CagdCrvFreeList(CubicCrvs);
  376.     return CubicCrvsMaxLen;
  377. }
  378.  
  379. /******************************************************************************
  380. * Given a bezier curve with order larger than cubic, approximate it using     *
  381. * piecewise cubic curves. A C^1 continuous approximation is computed so that  *
  382. * the approximation is at most sqrt(Tol2) away from the given curve.          *
  383. *   Input curve can be rational, although output is always a polynomial.      *
  384. ******************************************************************************/
  385. CagdCrvStruct *BzrApproxBzrCrvAsCubicPoly(CagdCrvStruct *Crv, CagdRType Tol2)
  386. {
  387.     CagdBType
  388.     ApproxOK = TRUE,
  389.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  390.     int i,
  391.     Order = Crv -> Order,
  392.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  393.     CagdPointType
  394.     PType = Crv -> PType,
  395.     CubicPType = CAGD_MAKE_PT_TYPE(FALSE, MaxCoord);
  396.     CagdCrvStruct *DistSqrCrv, *DiffCrv,
  397.     *CubicBzr = BzrCrvNew(4, CubicPType);
  398.     CagdRType E3Pt1[3], E3Pt2[3], E3Pt3[3], E3Pt4[3], Tan1[3], Tan2[3],
  399.     **CubicPoints = CubicBzr -> Points,
  400.     **Points = Crv -> Points;
  401.  
  402.     CagdCoerceToE3(E3Pt1, Points, 0, PType);
  403.     CagdCoerceToE3(E3Pt2, Points, 1, PType);
  404.     CagdCoerceToE3(E3Pt3, Points, Order - 2, PType);
  405.     CagdCoerceToE3(E3Pt4, Points, Order - 1, PType);
  406.     
  407.     /* End of the two points must match. */
  408.     for (i = 0; i < MaxCoord; i++) {
  409.     CubicPoints[i + 1][0] = E3Pt1[i];
  410.     CubicPoints[i + 1][3] = E3Pt4[i];
  411.     }
  412.  
  413.     /* Tangents at the end of the two points must match. */
  414.     PT_SUB(Tan1, E3Pt2, E3Pt1);
  415.     PT_SUB(Tan2, E3Pt4, E3Pt3);
  416.     PT_SCALE(Tan1, (Order - 1.0) / 3.0);
  417.     PT_SCALE(Tan2, (Order - 1.0) / 3.0);
  418.  
  419.     for (i = 0; i < MaxCoord; i++) {
  420.     CubicPoints[i + 1][1] = E3Pt1[i] + Tan1[i];
  421.     CubicPoints[i + 1][2] = E3Pt4[i] - Tan2[i];
  422.     }
  423.  
  424.     /* Compare the two curves by computed the distance square between */
  425.     /* corresponding parameter values.                      */
  426.     DiffCrv = CagdCrvSub(Crv, CubicBzr);
  427.     DistSqrCrv = CagdCrvDotProd(DiffCrv, DiffCrv);
  428.     CagdCrvFree(DiffCrv);
  429.     Points = DistSqrCrv -> Points;
  430.     if (IsNotRational) {
  431.     CagdRType
  432.         *Dist = Points[1];
  433.  
  434.     for (i = DistSqrCrv -> Order - 1; i >= 0; i--) {
  435.         if (*Dist++ > Tol2) {
  436.         ApproxOK = FALSE;
  437.         break;
  438.         }
  439.     }        
  440.     }
  441.     else {
  442.     CagdRType
  443.         *Dist0 = Points[0],
  444.         *Dist1 = Points[1];
  445.  
  446.     for (i = DistSqrCrv -> Order - 1; i >= 0; i--) {
  447.         if (*Dist1++ / *Dist0++ > Tol2) {
  448.         ApproxOK = FALSE;
  449.         break;
  450.         }
  451.     }        
  452.     }
  453.     CagdCrvFree(DistSqrCrv);
  454.  
  455.     if (ApproxOK)
  456.     return CubicBzr;
  457.     else {
  458.     CagdCrvStruct
  459.         *Crv1 = BzrCrvSubdivAtParam(Crv, 0.5),
  460.         *Crv2 = Crv1 -> Pnext,
  461.         *Crv1Approx = BzrApproxBzrCrvAsCubicPoly(Crv1, Tol2),
  462.         *Crv2Approx = BzrApproxBzrCrvAsCubicPoly(Crv2, Tol2);
  463.  
  464.     CagdCrvFree(Crv1);
  465.     CagdCrvFree(Crv2);
  466.  
  467.     for (Crv1 = Crv1Approx; Crv1 -> Pnext != NULL; Crv1 = Crv1 ->Pnext);
  468.     Crv1 -> Pnext = Crv2Approx;
  469.     return Crv1Approx;    
  470.     }
  471. }
  472.  
  473.